#include "./lj.h"

void mds_lj_init( mds_lj_data *lj, Real epsilon, Real sigma )
{
	assert(lj);
	memset( lj, 0, sizeof(mds_lj_data) );
	lj->sigma_sq     = sigma * sigma;
	lj->epsilon4     = MDS_F(4.0)  * epsilon;
	lj->epsilon24    = MDS_F(24.0) * epsilon;
}


Real mds_lj( const mds_vector *r, mds_vector *f, void *args )
{
	assert( r );
	assert( f );
	assert( args );
	{
		const mds_lj_data *lj     = (mds_lj_data *)args;
		const Real         r2     = mds_vector_normsq(r);
		const Real         q2     = lj->sigma_sq / r2; 
		const Real         q6     =  q2 * q2 * q2;
		const Real         q12    =  q6 * q6;

		/** compute force **/
		mds_vector_mulset(f, lj->epsilon24 * ( (q12+q12) - q6 ) / r2, r);
		
		/** return potential */
		return lj->epsilon4 * ( q12 -  q6 );
	}
}
